00100	comment
00200		this version of fake finds the angular and position
00300		errors for the edge operator as a function of step height,
00400		position, and angle;
00500	
00600	BEGIN "EDGE_TEST"
00700		comment
00800			find best choice of conf and diff for Hueckel edge
00900			operator based on chosen false positive rate for
01000			no signal.  For a given false positive rate, find
01100			choice which requires the least height step to
01200			achieve a 50% false negative rate.  This uses
01300			steps generated by a Monte Carlo procedure;
01400	
01500	DEFINE SAFE_ARRAY=1;	comment this makes all arrays safe;
01600	REQUIRE "VISION.2[SYS,TOB]" SOURCE_FILE; comment vision language package;
01700	REQUIRE EDGE_FILE;
01800	REQUIRE EDGE_PRINT;
01900	REQUIRE TV_VARIABLES;
02000	REQUIRE "STEP4.SAI[FAK,TOB]" SOURCE_FILE;
02100	
     

00100	comment
00200		conf=sqrt(s0sq/(s0sq+n0sq)) ie cosine
00300		diff=s0sq;
00400	
00500	BOOLEAN SIMPLE PROCEDURE YES(STRING S);
00600	BEGIN "YES"
00700	S←ASK_STRING(S);
00800	IF EQU(S,"Y")∨EQU(S,"y")∨EQU(S,"YES")∨EQU(S,"yes")
00900	THEN RETURN (1) ELSE RETURN (0);
01000	END "YES";
01100	DEFINE NBHD=8;	comment the size of the array with fake intensity values;
01200	LABEL FIRST_LOOP; ! entry for false positives and others;
01300	LABEL SECOND_LOOP; ! entry for false negatives;
01400	LABEL ERR_TEST; ! entry for angular and position errors;
01500	LABEL FALSE_NEG; ! entry for false negatives;
01600	LABEL AVG_DIFF;  ! entry for avg value of diff;
01700	LABEL HISTO_LOOP; ! loop entry for histograms repeat with change scale;
01800	DEFINE MAX_COUNT=2000; ! the max number of test samples generated;
01900	INTEGER_ARRAY INTENSITY[1:NBHD*NBHD]; ! store fake intensities here;
02000	REAL_ARRAY N0SQA[1:MAX_COUNT+20];	! store n0sq from edge oper here;
02100	REAL_ARRAY S0SQA[1:MAX_COUNT+20] ;	! store s0sq from edge oper here;
02200	INTEGER_ARRAY COUNT[0:100];	! histogram of scaled s0sq values;
02300	SHORT INTEGER K,K1;
02400	INTEGER	ENOUGH; ! how many samples (< max_count);
02500	SHORT REAL STEP_MEAN,STEP_AMPLITUDE,NOISE_STD;
02600		comment the mean,amplitude,noise of fake steps;
02700	SHORT REAL NSQ; ! scaling factor which depends on noise squared;
02800	SHORT REAL NOISE_SCALE;	! scaling factor for histograms;
02900	SHORT REAL THETA,RHO; ! angle,perp dist of edge from origin;
03000	SHORT INTEGER BYTES;	! number of bytes per word;
03100	SHORT REAL LAMBDA; ! wavelength (used to get phase near an edge);
03200	SHORT INTEGER N1,N5,KN1,KN5; comment 1% and 5% level counts and indices;
03300	
03400	ENOUGH←MAX_COUNT;
03500	
03600	comment now init the edge operator;
03700	CONF←.05;DIFF←.05; ! choose low values to reject nothing;
03800	BITS←6;
03900	MODIFY(BITS);
04000	LINLEN←WORDS(NBHD,BITS);
04100	TVWORD←GIOWD(INTENSITY);
04200	LSIDE←1;RSIDE←NBHD;
04300	FLINE←1;LLINE←NBHD;
04400	EJINIT;
04500	DIFF←0;CONF←.05;
04600	
04700	comment set up values for fake edges with zero height;
04800	NOISE_STD←5.0; ! amplitude of gaussian noise;
04900	STEP_MEAN←(2↑BITS-1)*.5; ! offset to get positive intensities;
05000	STEP_AMPLITUDE←0; ! step height for fake edges;
05100	NOISE_SCALE←20.0;
05200	
     

00100	comment find the angular and position errors at threshold;
00200	ERR_TEST:
00300	BEGIN "ERR_BLOCK"
00400	LABEL EDGE_APPLY;
00500	INTEGER K;
00600	INTEGER N;
00700	SHORT REAL STEP_FRACTION;
00800	SHORT REAL DIFF_FRACTION;
00900	SHORT REAL ALPHA,SUM_ALPHA,SUMSQ_ALPHA;
01000	SHORT REAL DELTA,SUM_DELTA,SUMSQ_DELTA;
01100	SHORT REAL STD_ALPHA,STD_DELTA;
01200	MODIFY(ENOUGH);
01300	MODIFY(STEP_FRACTION);
01400	MODIFY(NOISE_STD);
01500	STEP_AMPLITUDE←STEP_FRACTION*NOISE_STD;
01600	PRINT(STEP_AMPLITUDE);
01700	MODIFY(DIFF_FRACTION);
01800	DIFF←DIFF_FRACTION*NOISE_STD↑2;
01900	PRINT(DIFF);
02000	MODIFY(CONF);
02100	MODIFY(RHO);
02200	THETA←0;
02300	EDGE_APPLY:
02400		FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
02500		STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
02600			comment generate an array with zero edge,random;
02700		OP_APPLY(NBHD/2,NBHD/2); ! try the edge operator on it;
02800		ALPHA←ATAN2(CX,CY);
02900		TCRLF;PRINT(ALPHA);
03000		IF YES("TRY EDGE OPERATOR AGAIN? ") THEN GO TO EDGE_APPLY;
03100	SUM_ALPHA←SUMSQ_ALPHA←SUM_DELTA←SUMSQ_DELTA←0;
03200	N←0;
03300	FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ERR_LOOP"
03350	SHORT REAL THETA1;
03400	IF K MOD 20 = 0 THEN PRINT(K);
03450	THETA1←THETA+(π/4)*RANDOM-π/8; ! random angle +-π/8;
03500		FAKE(NBHD,NBHD,THETA1,RHO,LAMBDA,
03600		STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
03700			comment generate an array with zero edge,random;
03800		IF EJLI(NBHD/2,NBHD/2,0,0)∧¬BCOMP THEN BEGIN "INCR"
03900			N←N+1;
04000			ALPHA←ATAN2(CX,CY)-THETA1;
04100			SUM_ALPHA←SUM_ALPHA+ALPHA;
04200			SUMSQ_ALPHA←SUMSQ_ALPHA+ALPHA↑2;
04300			DELTA←SQRT((OPX-NBHD/2)↑2+(OPY-NBHD/2)↑2);
04400			SUM_DELTA←SUM_DELTA+DELTA;
04500			SUMSQ_DELTA←SUMSQ_DELTA+DELTA↑2;
04600			END "INCR";
04700	END "ERR_LOOP";
04800	TCRLF;PRINT(N);PRINT(SUMSQ_DELTA);PRINT(SUM_DELTA);
04900	TCRLF;PRINT(SUMSQ_ALPHA);PRINT(SUM_ALPHA);
05000	STD_ALPHA←SQRT((SUMSQ_ALPHA-SUM_ALPHA↑2/N)/N);
05100	STD_DELTA←SQRT((SUMSQ_DELTA-SUM_DELTA↑2/N)/N);
05200	TCRLF;
05300	PRINT(STD_ALPHA);
05400	PRINT(STD_DELTA);
05500	END "ERR_BLOCK";
05600	IF YES("REPEAT ERROR BLOCK? ") THEN GO TO ERR_TEST;
05700	END "EDGE_TEST";